Origins of agriculture density analysis

Current best version

Figure 1 - This figures shows…

Code

1. R environment setup
2. Setting time breaks
3. Defining origins
4. Import raster data
5. GAM smoothing models
6. Compare rates between origins and not-origins
7. Compare rates between different time periods
8. Setup final figure
9. Map for final figure
10. Trend through time panel for final figure
11. Assemble and print final figure

R environment setup

Attach libraries

library(png)
library(maptools)

Checking rgeos availability: TRUE

library(raster)
library(gam)

Loading required package: splines
Loading required package: foreach
foreach: simple, scalable parallel programming from Revolution Analytics
Use Revolution R for scalability, fault tolerance and more.
http://www.revolutionanalytics.com
Loaded gam 1.14

Set working directory

setwd("~/Desktop/Botero postdoc 2016/Human density and the origins of agriculture/")

Setting time breaks

Define the times of agricultural origins

par(mar=c(0,0,0,20))
d <- readPNG("Larson_dates.png")
plot(seq(0,18, length.out = 19), seq(0,36, length.out = 19), type="n",ylim=c(0,36),xlim=c(0, 18), xaxt="n")
rasterImage(d, 0,0,18,36, interpolate=TRUE, col=d)

Start_of_early_window <- 16-12
End_of_early_window_start_of_late_window <- 8.2
End_of_late_window <- 17-4.2
polygon(x=c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window), y=c(0, 34, 34, 0), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x=c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window), y=c(0, 34, 34, 0), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))

These dates are provided in the supplimentary information for the Larson (2014) paper. I’ve copied those values into a .csv table provided here.

domestication_times <- read.csv("Domestication timing larson 2014.csv")
dim(domestication_times)

[1] 77  8

Region	Species	Start.Exploitation	Finish.Exploitation	Start.predomestication	Finish.predomestication	Start.Domestication	Finish.Domestication
Southwest asia	Wheat	12.00	11.25	11.25	11.00	11.00	9.00
Southwest asia	Barley	12.00	11.25	11.25	10.50	10.50	9.00
Southwest asia	Lentil	12.00	11.00	11.00	10.50	10.50	9.00
Southwest asia	Pea	11.50	11.00	11.00	10.00	10.00	8.50
Southwest asia	Chickpea	11.00	10.50	10.50	10.25	10.25	8.25
Southwest asia	Broadbean	NA	NA	NA	NA	10.50	NA
Southwest asia	Flax	12.00	9.50	NA	NA	9.50	NA
Southwest asia	Olive	10.00	6.00	NA	NA	6.00	NA
Southwest asia	Sheep	12.00	10.50	10.50	9.75	9.75	8.00
Southwest asia	Goat	12.00	10.50	10.50	9.75	9.75	8.00
Southwest asia	Pig	12.00	11.50	11.50	9.75	10.25	9.00
Southwest asia	Cattle, taurine	11.50	10.50	10.50	10.25	10.25	8.00
Southwest asia	Cat	NA	NA	10.50	4.00	4.00	NA
South Asia	Tree cotton	8.50	4.50	NA	NA	4.50	NA
South Asia	Rice	8.00	5.00	5.00	4.00	4.00	2.50
South Asia	Little millet	NA	NA	NA	NA	4.50	NA
South Asia	Browntop millet	NA	NA	NA	NA	4.00	NA
South Asia	Mungbean	NA	NA	4.50	3.50	3.50	3.00
South Asia	Pigeonpea	NA	NA	NA	NA	3.50	NA
South Asia	Zebu cattle	9.00	8.00	NA	NA	8.00	6.50
South Asia	Water buffalo	6.00	4.50	NA	NA	4.50	NA
East Asia	Broomcorn Millet	10.00	8.00	NA	NA	8.00	NA
East Asia	Foxtail millet	11.50	7.50	NA	NA	7.50	NA
East Asia	Rice	10.00	8.00	8.00	7.50	7.50	5.00
East Asia	Soybean	8.50	5.50	NA	NA	5.50	4.00
East Asia	Ramie	NA	NA	NA	NA	5.25	NA
East Asia	Melon	7.00	4.00	NA	NA	4.00	3.75
East Asia	Pig	12.00	8.50	NA	NA	8.50	6.00
East Asia	Silkworm	7.00	5.25	NA	NA	5.25	NA
East Asia	Yak	NA	NA	NA	NA	4.25	NA
East Asia	Horse	7.50	6.75	6.75	5.50	5.50	4.00
East Asia	Bactrian Camel	NA	NA	NA	NA	4.50	NA
East Asia	Duck	2.50	1.00	NA	NA	1.00	NA
East Asia	Chicken	6.00	4.00	NA	NA	4.00	NA
New Guinea	Banana	10.00	7.00	7.00	4.00	4.00	NA
New Guinea	Taro	10.00	7.00	7.00	4.00	NA	NA
New Guinea	Yam	10.00	7.00	7.00	4.00	NA	NA
Africa and Arabia	Date palm	7.00	6.00	NA	NA	5.00	NA
Africa and Arabia	Sorghum	8.00	4.00	NA	NA	4.00	NA
Africa and Arabia	Pearl millet	NA	NA	NA	NA	4.50	3.50
Africa and Arabia	Fonio	NA	NA	NA	NA	2.50	NA
Africa and Arabia	Cowpea	NA	NA	NA	NA	3.75	NA
Africa and Arabia	Hyacinth bean	NA	NA	NA	NA	3.75	NA
Africa and Arabia	Rice	3.50	2.00	NA	NA	2.00	NA
Africa and Arabia	Oil palm	9.25	3.50	NA	NA	3.50	NA
Africa and Arabia	Cattle	NA	NA	9.00	7.75	7.75	6.50
Africa and Arabia	Donkey	9.00	5.50	NA	NA	5.50	3.50
Africa and Arabia	Dromedary camel	6.50	3.00	NA	NA	3.00	NA
Africa and Arabia	Guinea fowl	NA	NA	2.50	1.50	1.50	NA
North America	Squash	6.50	5.00	NA	NA	5.00	NA
North America	Sunflower	6.00	4.75	NA	NA	4.00	NA
North America	Sumpweed	6.00	4.50	NA	NA	4.00	NA
North America	Pitseed goosefoot	4.75	3.75	NA	NA	3.75	NA
Meso-america	Squash (pepo)	NA	NA	NA	NA	10.00	9.50
Meso-america	Maize	10.00	9.00	NA	NA	9.00	NA
Meso-america	Foxtail millet-grass	NA	NA	NA	NA	6.00	4.00
Meso-america	Common bean	NA	NA	NA	NA	3.00	NA
Meso-america	Avocado	NA	NA	NA	NA	3.00	NA
Meso-america	Chile pepper	NA	NA	NA	NA	3.00	NA
Meso-america	Turkey	NA	NA	NA	NA	2.00	NA
South America	Chili pepper	NA	NA	NA	NA	6.00	NA
South America	Peanut	NA	NA	NA	NA	5.00	NA
South America	Cotton	NA	NA	NA	NA	6.00	NA
South America	Coca	NA	NA	NA	NA	8.00	NA
South America	Now-minor root crops (arrowroot, leren)	NA	NA	NA	NA	9.00	NA
South America	Squash	NA	NA	NA	NA	10.00	NA
South America	Common bean	NA	NA	NA	NA	5.00	NA
South America	Lima bean	NA	NA	8.25	NA	6.00	NA
South America	Monioc	NA	NA	NA	NA	7.00	NA
South America	Sweet potato	NA	NA	NA	NA	5.00	NA
South America	White potato	7.00	4.50	NA	NA	4.00	NA
South America	Quinoa	5.00	NA	NA	NA	3.50	NA
South America	Yam	NA	NA	NA	NA	5.50	NA
South America	Llama	10.00	6.00	NA	NA	6.00	NA
South America	Alpaca	10.00	5.00	NA	NA	5.00	NA
South America	Guinea pig	NA	NA	NA	NA	5.00	NA
South America	Muscovy Duck	NA	NA	NA	NA	4.00	NA

par(mar=c(5,4,6,1))
dates <- unlist(domestication_times[3:8])
hist(dates, breaks = 22, xlim=c(15,0), xlab="K years ago", col=adjustcolor("cornflowerblue", alpha= 0.5), border=adjustcolor("cornflowerblue", alpha= 0.9), main="All dates in dataset"  )
mtext("This tells us about how evenly our evidence is distributed in time", 3, line=1)

hist(dates, breaks = 22, xlim=c(15,0), xlab="Thousand years ago", col=adjustcolor("cornflowerblue", alpha= 0.5), border=adjustcolor("cornflowerblue", alpha= 0.9), main="All dates in dataset with Larson(2014) date windows")
Start_of_early_window <- 12
End_of_early_window_start_of_late_window <- 8.2
End_of_late_window <- 4.2
polygon(x=c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window), y=c(0, 30, 30, 0), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))

polygon(x=c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window), y=c(0, 30, 30, 0), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
hist(dates, breaks = 22, xlim=c(15,0), xlab="K years ago", col=adjustcolor("cornflowerblue", alpha= 0.2), border=adjustcolor("cornflowerblue", alpha= 0.9), main="", add=TRUE)

mtext("Early Holocene", 3, line = -1, adj=.3)
mtext("Middle Holocene", 3, line= -1, adj=.6)

par(mfrow=c(2,3), mar=c(4,4,2,0))
dim(domestication_times)

[1] 77  8

specific_dates <- domestication_times[,3:8]
for(i in c(1, 3, 5, 2, 4, 6)){
hist(specific_dates[,i], breaks = 22, xlim=c(15,0), xlab="Thousand years ago", col=adjustcolor("cornflowerblue", alpha= 0.5), border=adjustcolor("cornflowerblue", alpha= 0.9), main= names(specific_dates)[i])
Start_of_early_window <- 12
End_of_early_window_start_of_late_window <- 8.2
End_of_late_window <- 4.2
polygon(x=c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window), y=c(0, 30, 30, 0), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x=c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window), y=c(0, 30, 30, 0), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
hist(specific_dates[,i], breaks = 22, xlim=c(15,0), xlab="K years ago", col=adjustcolor("cornflowerblue", alpha= 0.2), border=adjustcolor("cornflowerblue", alpha= 0.9), main="", add=TRUE)
}

I’m creating new rows for this table, combining dates in different ways to make the CDFs below look more authentic. This makes it so that pre-ag always happens before post-ag. What I’ve done is given the later date to the earlier date when those dates are missing.

h <- which(is.na(domestication_times[,3]))
domestication_times <- cbind(domestication_times, rep(NA, length(domestication_times[,1])))
domestication_times[,9] <- domestication_times[,3]
domestication_times[h,9] <- domestication_times[h,7]
colnames(domestication_times)[9] <- "adopt exploitation date"
domestication_times[,10] <- domestication_times[,7]
domestication_times[which(is.na(domestication_times[,10])),10] <- 0
colnames(domestication_times)[10] <- "start of ag"
#save(domestication_times, file="~/Desktop/Human density and the origins of agriculture/Domestication timing larson 2014.Rdata")

I think these are best described by a cummulative distribution, showing how they accumulate over time.

for(i in 1:8){
type_number <- i
    match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 9]
    maxer <- max(match, na.rm=TRUE)
    j <- ecdf(maxer-match)
    print(levels(domestication_times$Region)[ type_number])
    print(match)
    print(j)
}

[1] "Africa and Arabia"
 [1] 7.00 8.00 4.50 2.50 3.75 3.75 3.50 9.25 7.75 9.00 6.50 1.50
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:11] =      0,   0.25,   1.25,  ...,   6.75,   7.75
[1] "East Asia"
 [1] 10.00 11.50 10.00  8.50  5.25  7.00 12.00  7.00  4.25  7.50  4.50  2.50  6.00
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:11] =      0,    0.5,      2,  ...,   7.75,    9.5
[1] "Meso-america"
[1] 10 10  6  3  3  3  2
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:4] =      0,      4,      7,      8
[1] "New Guinea"
[1] 10 10 10
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:1] =      0
[1] "North America"
[1] 6.50 6.00 6.00 4.75
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:3] =      0,    0.5,   1.75
[1] "South America"
 [1]  6.0  5.0  6.0  8.0  9.0 10.0  5.0  6.0  7.0  5.0  7.0  5.0  5.5 10.0 10.0  5.0  4.0
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:8] =      0,      1,      2,  ...,      5,      6
[1] "South Asia"
[1] 8.5 8.0 4.5 4.0 3.5 3.5 9.0 6.0
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:7] =      0,    0.5,      1,  ...,      5,    5.5
[1] "Southwest asia"
 [1] 12.0 12.0 12.0 11.5 11.0 10.5 12.0 10.0 12.0 12.0 12.0 11.5  4.0
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:6] =      0,    0.5,      1,  ...,      2,      8

par(mfcol=c(2,5), mar=c(4,0,5,0))
plot(0,0, type="n", xaxt="n", xlab="", bty="n")
mtext("Percent of species that will eventually \n be domesticated in a region", 2, line=-5, cex=0.5)

plot(0,0, type="n", xaxt="n", xlab="", bty="n")
mtext("Percent of species that will eventually \n be domesticated in a region", 2, line=-5, cex=0.5)

for(i in 1:8){
type_number <- i
    match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 9]
    maxer <- max(match, na.rm=TRUE)
    j <- ecdf(maxer-match)
    #print(j)
    
plot(0,0, xlim=c(15,0), ylim=c(0,100), ylab="Percent of species that will eventually \n be domesticated in a region", xlab="Thousand years ago", main=levels(domestication_times$Region)[ type_number], type="n", yaxt="n")
x_seq <- rev(c(0,seq(0, maxer, length.out=100)))
y_seq <- 100 * (c(0, j(seq(0, maxer, length.out=100))))
lines(x_seq, y_seq,  ylim=c(-1,1))
polygon(c(0, x_seq), c(0, y_seq), border=adjustcolor("cornflowerblue",alpha=1), col=adjustcolor("cornflowerblue", alpha=0.2))
if(i == 2 | i == 1)axis(2)
if(i == 3)mtext("Cummulative distribution function for the accumulation of domesticates", 3, line=3.8, col="cornflowerblue")
}

par(mfcol=c(2,5), mar=c(4,0,5,0))
plot(0,0, type="n", xaxt="n", xlab="", bty="n")
mtext("Percent of species that will eventually \n be domesticated in a region", 2, line=-5, cex=0.5)

plot(0,0, type="n", xaxt="n", xlab="", bty="n")
mtext("Percent of species that will eventually \n be domesticated in a region", 2, line=-5, cex=0.5)

for(i in 1:8){
type_number <- i
    match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 9]
    maxer <- max(match, na.rm=TRUE)
    j <- ecdf(maxer-match)
    #print(j)
    
plot(0,0, xlim=c(15,0), ylim=c(0,100), ylab="Percent of species that will eventually \n be domesticated in a region", xlab="Thousand years ago", main=levels(domestication_times$Region)[ type_number], type="n", yaxt="n")
x_seq <- rev(c(0,seq(0, maxer, length.out=100)))
y_seq <- 100 * (c(0, j(seq(0, maxer, length.out=100))))
lines(x_seq, y_seq,  ylim=c(-1,1))
polygon(c(0, x_seq), c(0, y_seq), border=adjustcolor("cornflowerblue",alpha=1), col=adjustcolor("cornflowerblue", alpha=0.2))
abline(v= maxer - quantile(j)[2], col="limegreen", lwd=2)
if(i == 2 | i == 1)axis(2)
if(i == 2)mtext("25%", 3, line=3.5, adj=-1, col="limegreen")
if(i == 3)mtext("Cummulative distribution function for the accumulation of domesticates", 3, line=3.8, col="cornflowerblue")
if(i == 4)mtext("Choose a y to predict an x", 3, line=3.3, col="cornflowerblue")
    break_one <- maxer
            break_two <- maxer - quantile(j)[2]
                
    polygon(x=c(break_two, break_two, break_one, break_one), y=c(0, 1, 1, 0), col=adjustcolor("cornflowerblue", alpha=0.2), border=adjustcolor("cornflowerblue",alpha=1))
            lines(x=c(break_two, break_two), y=c(0,-1), col="cornflowerblue")
            abline(h = 25, col="limegreen", lwd=2)
}

Make this a function. There is a choice of two methods here. At the end of this section we need to print the desision we’re passing to the later analyses.

Defining origins

origins <- readShapePoly('Origins_updated.shp')
origin.time.region <- c(2, 2, 1, 1, 1, 2, 2, 1, 2, 2, 
                        2, 2, 1, 2, 2, 2, 2, 2, 2, 2) # 1 = early; 2 = middle
as.character(origins$CONTINENT)

 [1] "W_African_Sav" "Sudanic_Savan" "West Africa T" "Ethipian plat" NA             
 [6] "E_North_Ameri" "New_Guinea"    "Mesoamerica"   "N_Lowland_SA"  "NW_Lowland_SA"
[11] "Sava_W_India"  "S_India"       "Ganges_E_Indi" "Chinese_loess" "Japanese"     
[16] "Lower-MiddleY" "South trop ch" NA              "Southwes amaz" "C/S_Andes"    

#subset_order <- c(1, 2, 3, 5, 6, 8, 9, 10, 11, 12, 17, 18)
subset_order <- c(8, 10, 9, 5, 18, 7, 6, 20, 1, 2, 13, 16)
origins_subset <- origins[subset_order,]
origins_subset$CONTINENT

 [1] Mesoamerica   NW_Lowland_SA N_Lowland_SA  <NA>          <NA>          New_Guinea   
 [7] E_North_Ameri C/S_Andes     W_African_Sav Sudanic_Savan Ganges_E_Indi Lower-MiddleY
18 Levels: C/S_Andes Chinese_loess E_North_Ameri Ethipian plat Ganges_E_Indi ... West Africa T

origins_subset$name

NULL

library(maps)
map()
map(origins, add=TRUE, fill=TRUE, col=adjustcolor("cornflowerblue", alpha=1))

database does not (uniquely) contain the field 'name'.

map()
d <- readPNG("Larson_origins.png")
rasterImage(d, -180, -90, 180, 110, interpolate=TRUE, col=d)

map(add=TRUE)
map(origins, add=TRUE, fill=TRUE, col=adjustcolor("cornflowerblue", alpha=1))

database does not (uniquely) contain the field 'name'.

# need to reproject

This is obviously a bad projection fit right now.

Import raster data

#subset and reorder origins. This is currently done at the end of the plot but should be moved forward.
# Load data for population density
load("PopD_all_December.rdata")
PopD.ALL

class       : RasterStack 
dimensions  : 288, 720, 207360, 18  (nrow, ncol, ncell, nlayers)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -60, 84  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
names       :        fourK,        fiveK,         sixK,       sevenK,       eightK,        nineK,         tenK,      elevenK,      twelveK,    thirteenK,    fourteenK,     fifteenK,     sixteenK,   seventeenK,    eighteenK, ... 
min values  : 5.611358e-07, 1.067142e-06, 2.508241e-06, 6.317553e-06, 2.286934e-05, 7.631922e-05, 1.272693e-04, 2.118215e-04, 2.602175e-04, 3.226203e-04, 4.390267e-04, 5.572032e-04, 7.313966e-04, 8.286005e-04, 8.297062e-04, ... 
max values  :     2.051069,     2.013452,     2.142908,     1.888403,     1.863014,     1.880628,     1.650615,     1.678033,     1.697732,     1.499115,     1.517264,     1.443677,     1.464867,     1.453581,     1.436394, ... 

# Extract data to a matrix
Pop <- values(PopD.ALL)
r <- raster(PopD.ALL, 1)
r

class       : RasterLayer 
dimensions  : 288, 720, 207360  (nrow, ncol, ncell)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -60, 84  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : fourK 
values      : 5.611358e-07, 2.051069  (min, max)

GAM smoothing models

Justification for General Adative Models.

We need to justify our decision to use a GAM over other models. This should include citations to back up those arguments.

Fit and plot GAM model with different degrees of freedom

We should make our decisions very transparent here. We should be able to justify our decision of 3 degrees of freedom over other possible values.

Density projections

# Read the polygons
origins <- readShapePoly('Origins_updated.shp')

# Extract data
per.origin <- extract(r, origins, cellnumber = TRUE, buffer = 100000)

NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files

names(per.origin) <- origins@data[, 1]
str(per.origin)

List of 20
 $ W_African_Sav: num [1:309, 1:2] 92506 92507 92508 92509 92510 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:2] "cell" "value"
 $ Sudanic_Savan: num [1:306, 1:2] 99050 99051 99052 99053 99054 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:2] "cell" "value"
 $ West Africa T: num [1:427, 1:2] 102609 102610 102611 103326 103327 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:2] "cell" "value"
 $ Ethipian plat: num [1:275, 1:2] 106281 106282 106283 106998 106999 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:2] "cell" "value"
 $ NA           : num [1:195, 1:2] 67404 67405 67406 68119 68120 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:2] "cell" "value"
 $ E_North_Ameri: num [1:170, 1:2] 64270 64271 64984 64985 64986 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:2] "cell" "value"
 $ New_Guinea   : num [1:24, 1:2] 127360 127361 127362 128080 128081 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:2] "cell" "value"
 $ Mesoamerica  : num [1:73, 1:2] 93032 93033 93034 93035 93036 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:2] "cell" "value"
 $ N_Lowland_SA : num [1:39, 1:2] 110373 110374 110375 111092 111093 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:2] "cell" "value"
 $ NW_Lowland_SA: num [1:24, 1:2] 123319 123320 123321 123322 123323 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:2] "cell" "value"
 $ Sava_W_India : num [1:34, 1:2] 83312 84031 84032 84749 84750 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:2] "cell" "value"
 $ S_India      : num [1:18, 1:2] 99153 99154 99872 99873 99874 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:2] "cell" "value"
 $ Ganges_E_Indi: num [1:92, 1:2] 85494 85495 85496 85497 85498 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:2] "cell" "value"
 $ Chinese_loess: num [1:84, 1:2] 72598 72599 73318 73319 73320 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:2] "cell" "value"
 $ Japanese     : num [1:36, 1:2] 59681 59682 60401 61121 61843 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:2] "cell" "value"
 $ Lower-MiddleY: num [1:131, 1:2] 72547 72548 72549 73267 73268 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:2] "cell" "value"
 $ South trop ch: num [1:178, 1:2] 84818 84819 84820 84821 84822 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:2] "cell" "value"
 $ NA           : num [1:258, 1:2] 62493 62494 62495 62496 62497 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:2] "cell" "value"
 $ Southwes amaz: num [1:194, 1:2] 137758 137759 137760 137761 137762 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:2] "cell" "value"
 $ C/S_Andes    : num [1:165, 1:2] 137736 137737 137738 137739 137740 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:2] "cell" "value"

par(mfrow=c(4,5), mar=c(0,0,0,0))
for(h in 1:20){
#h <- 3
originI <- Pop[per.origin[[h]][, 1], ]
x_values <- matrix(c(4:21), dim(originI)[1], 18, byrow=TRUE)
x_value_vector <- as.vector(x_values)
y_value_vector <- scale(as.vector(originI))
density_trend <- cbind(x_value_vector, y_value_vector)
plot(density_trend, col=adjustcolor("cornflowerblue", alpha=0.8),  xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(-10,10), xlim=c(21,4))
points(density_trend, cex=0.5, col=adjustcolor("cornflowerblue", 0.5))
}

par(mfrow=c(4,5), mar=c(0,0,0,0))
for(h in 1:20){
#h <- 3
originI <- Pop[per.origin[[h]][, 1], ]
x_values <- matrix(c(4:21), dim(originI)[1], 18, byrow=TRUE)
x_value_vector <- as.vector(x_values)
y_value_vector <- scale(as.vector(originI))
density_trend <- cbind(x_value_vector, y_value_vector)
plot(density_trend, col=adjustcolor("cornflowerblue", alpha=0.8),  xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(-4,8), xlim=c(21,4))
points(density_trend, cex=0.5, col=adjustcolor("grey", 0.1))
gammer <- loess(density_trend[,2] ~ density_trend[,1] )
summary(gammer)
lines(density_trend[,1] ,predict(gammer),  col="cornflowerblue", lwd=2)
}

span
the parameter α which controls the degree of smoothing.

The size of the neighbourhood is controlled by α (set by span or enp.target). For α < 1, the neighbourhood includes proportion α of the points, and these have tricubic weighting (proportional to (1 - (dist/maxdist)³⁾3). For α > 1, all points are used, with the ‘maximum distance’ assumed to be α^(1/p) times the actual maximum distance for p explanatory variables.

par(mfrow=c(4,5), mar=c(0,0,0,0))
for(h in 1:20){
#h <- 3
originI <- Pop[per.origin[[h]][, 1], ]
x_values <- matrix(c(4:21), dim(originI)[1], 18, byrow=TRUE)
x_value_vector <- as.vector(x_values)
y_value_vector <- scale(as.vector(originI))
density_trend <- cbind(x_value_vector, y_value_vector)
#points(density_trend, cex=0.5)
gammer <- loess(density_trend[,2] ~ density_trend[,1], span = .5)
summary(gammer)
predict_gam <- predict(gammer, se=TRUE)
plot(density_trend, col=adjustcolor("cornflowerblue", alpha=0.8),  xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(-1.35,1.35), xlim=c(21,4))
polygon(x=c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window), y=c(-2, 2, 2, -2), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x=c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window), y=c(-2, 2, 2, -2), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
abline(h=0, lty=2, col="grey")
lines(density_trend[,1] ,predict_gam$fit,  col="cornflowerblue")
lines(density_trend[,1] ,  predict_gam$fit + predict_gam$se.fit,  col="grey", lty=1)
lines(density_trend[,1] ,  predict_gam$fit - predict_gam$se.fit,  col="grey", lty=1)
}

# need to add a global mean, an everything but the origins mean, and a buffer around the origins mean. 
# Function standardization
std <- function(x) {
  b <- (x - min(x)) / (max(x) - min(x))
  return(rev(b))
}
diff_df <- function(h){ 
# Calculating mean and 
global.means <- global.SD <- list()
for (j in 1:length(per.origin)) {
  #print(j)
  originI <- Pop[per.origin[[j]][, 1], ]
  time <- 21:4
  originI <- na.exclude(originI)
  b <- apply(originI, 1, std)
  nJ <- nrow(originI)
  predictions <- matrix(nrow = nJ, ncol = length(time))
  colnames(predictions) <- as.character(time)
  for(i in 1:nJ) {
    
    # Need to show a gradient of these df values. 
    model <- gam(b[, i] ~ s(time, df = h))
    col <- sample(rainbow(100), 1)
    predictions[i, ] <- predict(model)
  }
  global.means[[j]] <- apply(predictions, 2, mean) 
  global.SD[[j]] <- apply(predictions, 2, sd)
}
names(global.means) <- paste(names(per.origin), "Means")
names(global.SD) <- paste(names(per.origin), "SD")
return(list(global.means, global.SD))
}

# Confirm that the x-axis is oriented correctly
for_3 <- diff_df(1)
global.means <- for_3[[1]]
global.SD <- for_3[[2]]
i <- 1
plot(4:21, global.means[[i]], col=adjustcolor("cornflowerblue", alpha=0.8),  xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(0,1), xlim=c(0,20))
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))

polygon(x= 22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
originI <- Pop[per.origin[[1]][, 1], ]
#colnames(originI)
#plot(21:4, rep(0, length(21:4)), type="n", ylim=c(0,1), xlim=c(21,4), xaxt="n", yaxt="n")
for(j in 1:18){
points(23- c(jitter(rep(4, length(originI[,j])), 4.5)  + j), originI[,j], pch=19, cex=.1)
}

polygon(c(4:21,21:4) -3 , c(global.means[[i]] + abs(global.SD[[i]]), rev(global.means[[i]] - abs(global.SD[[i]]))), col="cornflowerblue")
lines(global.means[[i]], type="b", pch=names(global.means[[i]]))

axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
mtext(names(global.means)[1], 3, line=-2)

means_matrix <- matrix(rep(NA,19*20), 20, 19)
colnames(means_matrix) <- c("origin", rev(seq(4, 21, by=1)))
means_matrix[,1] <- names(global.means)
for(i in 1:20){
means_matrix[i,2:19] <- global.means[[i]]
}

kable(means_matrix, caption= "Mean values")

origin	21	20	19	18	17	16	15	14	13	12	11	10	9	8	7	6	5	4
W_African_Sav Means	0.275750492966357	0.237971703906143	0.200504339081811	0.165047402018049	0.135516039740239	0.121748343406075	0.132429431895417	0.17103635569842	0.235948736377181	0.321617048802923	0.418964971526823	0.514676226549989	0.591121916502731	0.628378724452797	0.621163452983342	0.582047087468811	0.532681459452571	0.482909221409947
Sudanic_Savan Means	-0.0186911478417769	-0.00828615401443568	0.00413446846041128	0.0215161762523414	0.0473073567212149	0.0847812201109172	0.136209309975754	0.20228390605606	0.282347797217804	0.374044632344897	0.472296827728429	0.568020597195695	0.648225465565203	0.697539021870325	0.713328901689326	0.707327373243822	0.697887725176853	0.69056153763021
West Africa T Means	-0.0145401490853128	-0.00333603963380315	0.00977013345024412	0.027586956435481	0.0533165943413057	0.0895594750479319	0.138174778488521	0.199879802842878	0.274517388507917	0.360833593695506	0.45568720488448	0.552897898880818	0.643252356019774	0.715741781485123	0.768806992435192	0.809952043255802	0.850519193881187	0.893430170662931
Ethipian plat Means	0.0516408507031156	0.0622592531360626	0.0750787355293174	0.0933555709623621	0.12092332890413	0.16140265228906	0.217329358000213	0.28895304336974	0.373944304902601	0.467450912361356	0.562582640478741	0.650493007468181	0.722074582042064	0.768618080463592	0.790630347643354	0.797506826970228	0.803240018411344	0.812789660231236
NA Means	0.72808776413489	0.6723269776978	0.613165929035697	0.547108730705425	0.472515896153583	0.393590592158964	0.319669766181732	0.262708374009385	0.233548745835161	0.239616095881227	0.28323109870948	0.359870592475047	0.454992276384195	0.545528134731924	0.616545768777106	0.664148269449493	0.691583655592282	0.708950521834131
E_North_Ameri Means	-0.0424144312128647	-0.0202618739720356	0.00475194649595462	0.0365915191653316	0.079541401163135	0.136467523436837	0.207873937132646	0.291025513809716	0.380397583621781	0.468564149521066	0.547361392244133	0.609663875755146	0.654835785797328	0.689283124023394	0.721880145858871	0.761679332218606	0.814909663949399	0.875153810181115
New_Guinea Means	0.304105470512681	0.28132745701113	0.25734789630324	0.230996250201583	0.201793192056862	0.171544911503563	0.144548456661618	0.127184052406642	0.126877035836262	0.15052865629083	0.203062414319868	0.285273581529037	0.389014353820968	0.496874952215665	0.597210277530265	0.683771044637899	0.754313180721177	0.816774228406938
Mesoamerica Means	-0.0315889824677606	-0.00787599693194719	0.0184956204628149	0.0511692835844435	0.0938250236223921	0.148344738459587	0.214601573425648	0.290380258066549	0.372334867390379	0.456479966860125	0.538354702268212	0.613001244481024	0.676612820283987	0.726625668176205	0.767020170677251	0.80679724213028	0.855043623685109	0.908351302981932
N_Lowland_SA Means	0.208792535197897	0.197230324856633	0.185079121786267	0.171696272010441	0.156804931838545	0.141342946965234	0.128176787393512	0.122340507893562	0.130435008736839	0.159548596595071	0.215485053059509	0.300036154887671	0.405861213383308	0.51607541749585	0.62000422933541	0.712583747106496	0.792272847070064	0.866211062729842
NW_Lowland_SA Means	0.326576386801442	0.312436135756889	0.296990439986673	0.278723633997093	0.256416348833303	0.230169366346516	0.201945873000531	0.175714919736147	0.157651360929879	0.155352073191265	0.176275702224548	0.225642677767534	0.302153168560349	0.396329228756827	0.497310353388214	0.594432476311419	0.67925446713634	0.756553171147025
Sava_W_India Means	-0.0237936944828927	-0.00818464775151594	0.00933560457136895	0.0313826871867378	0.0607897006722834	0.0992817835941126	0.147568966381951	0.205468601119744	0.272443196878899	0.347413657258897	0.428584371379107	0.513039004997093	0.595722652725526	0.669265301648986	0.731180489316986	0.788223314024364	0.850666362502511	0.917274849593167
S_India Means	-0.0222028081854339	-0.0066002592031116	0.0109240264309662	0.0332350046742644	0.0636053304857559	0.104248025366543	0.156161312176903	0.218860716380081	0.291071667647164	0.371586637639489	0.458699358753368	0.549398894270441	0.638133144838091	0.71736832702396	0.786036841922908	0.849257343160103	0.913846434565112	0.979674443478567
Ganges_E_Indi Means	-0.0168749067644347	-0.00665207444708958	0.00516238304586039	0.0208799296723999	0.0431774421109922	0.0742138287818073	0.115649991879242	0.168563999805984	0.233555951882637	0.310487513335769	0.398078302603856	0.492768744130255	0.588625507090279	0.676929905267031	0.753772199676962	0.823106046514196	0.892776064568763	0.964026386513837
Chinese_loess Means	-0.00992922030088753	-0.00685627234272808	-0.00263644972242058	0.00456650135793812	0.0171800768542588	0.0379487356429578	0.0695048935574124	0.114680636092733	0.176234630227922	0.255681643811977	0.353044690591995	0.465903186889885	0.587480671879438	0.70545496017147	0.811956397518547	0.90571144015868	0.990020919073498	1.06943656618837
Japanese Means	0.474649163724578	0.443432703865345	0.410041741204845	0.372273786700216	0.329121724356353	0.283030955690219	0.240090867441392	0.209371090119366	0.201589325186606	0.226840139404003	0.291350184139259	0.39377817553382	0.519726499901562	0.642890482156585	0.747382040461932	0.829407324688329	0.892195580878682	0.945280306169006
Lower-MiddleY Means	-0.00682865849885968	-0.00060483620796883	0.00671174779516591	0.0168353119264841	0.0318371795131959	0.0536523053498361	0.0841840911894777	0.125268923860985	0.179000307645022	0.247353016268821	0.331229454779227	0.42919493710873	0.535517759646172	0.639543869787761	0.736041145928061	0.82613958614389	0.915048715480776	1.00384399368775
South trop ch Means	-0.00259675023606447	-0.00369446566682659	-0.00363747351558766	-0.000390054740314468	0.00904019752001561	0.0284916546492706	0.0620623767539311	0.113139640613902	0.183967610626003	0.275388397726442	0.385135385537217	0.507968233105058	0.634326208132707	0.750251125909705	0.848539337462573	0.930537024584384	1.00231564635787	1.06963296926003
NA Means	0.0366402042172722	0.0345449445801444	0.0333093635431334	0.0345629129341452	0.0409323268055525	0.0561025771005159	0.0840836906480487	0.12805515000113	0.189442050381472	0.269535827938253	0.368930130834347	0.484509427569726	0.606456684388843	0.720252238756432	0.817753922638379	0.899780585526795	0.971748613294788	1.03906371492767
Southwes amaz Means	0.207835723880993	0.197019933600991	0.185689737006914	0.173364610555289	0.159843367441753	0.145872982114604	0.134154420414965	0.129609533908725	0.138844389475622	0.168981654647338	0.225948263430072	0.312180448593161	0.422378015176177	0.542599962121189	0.664662002239541	0.784132607465062	0.897503322879635	1.00679099727573
C/S_Andes Means	0.0291205853249027	0.0314816768181023	0.0344713731256695	0.0390464429788137	0.0464560633526933	0.0581169360685076	0.0760540915093959	0.103031413769306	0.142906851931112	0.200040121010661	0.27794292932455	0.377058478195417	0.491456226460239	0.608384004959301	0.721451798462705	0.828972489197	0.930712639008333	1.02909296975965

#global.SD
SD_matrix <- matrix(rep(NA,19*20), 20, 19)
colnames(SD_matrix) <- c("origin", rev(seq(4, 21, by=1)))
SD_matrix[,1] <- names(global.SD)
for(i in 1:20){
SD_matrix[i,2:19] <- global.SD[[i]]
}

kable(SD_matrix, caption= "SD values")

origin	21	20	19	18	17	16	15	14	13	12	11	10	9	8	7	6	5	4
W_African_Sav SD	0.186556185173684	0.159071135711959	0.130974576751884	0.101946183025733	0.0724951714757066	0.0464817364491077	0.0281882850407625	0.0212718909296576	0.0230484233864053	0.0260000666669227	0.0272778672869396	0.0272414655989988	0.0286145237100166	0.0359445310006679	0.0509667284519141	0.0718957662099432	0.0961865477719618	0.12156398251189
Sudanic_Savan SD	0.0182791317887036	0.0120186235797875	0.00857406807198337	0.0119997090375925	0.0200097740504662	0.0296208384871706	0.0390634218166876	0.0465027632519908	0.0502676895837266	0.0490728899403298	0.0422563308391379	0.0301155921728213	0.016251502147543	0.0230519738470708	0.0490762671854076	0.0798060762063573	0.111190090221826	0.142489654117473
West Africa T SD	0.0128551006323539	0.0100809743903953	0.0110855392306364	0.0157356177847031	0.0225339267643659	0.0302798766216314	0.0378753478987503	0.0441002370122807	0.0478360621585905	0.0482472418609378	0.045098717097735	0.0390673319501593	0.0318512113726889	0.0254461299782974	0.0263553160311764	0.0393349506349416	0.0572833097318681	0.0763495674818114
Ethipian plat SD	0.198659638310769	0.170440069868076	0.141737372134567	0.112259279537512	0.0864755397538964	0.0791421681652602	0.0986717402217165	0.129156174657549	0.154623472712395	0.166247254796878	0.159947337518384	0.135808529359533	0.0994473010536208	0.0618393023227622	0.0356919120363483	0.0377802530877869	0.0546872860266305	0.0726120006791913
NA SD	0.352588159054744	0.332443095173091	0.310979361000686	0.28580375866142	0.254547237522766	0.216346143796596	0.171980702025272	0.12388521679558	0.0763676236214709	0.0403854093134876	0.0435516738756018	0.0673123480707615	0.0842474074682211	0.089847619869054	0.0865158533690929	0.0824620916709309	0.0858640772961775	0.09604764429718
E_North_Ameri SD	0.00744774460623942	0.00688305183133896	0.00730077417171542	0.00859083907761313	0.0105320687647222	0.012800810910875	0.0151914913730156	0.0175920880323135	0.0199689872021909	0.0222421032252577	0.0241465903979098	0.0252306084926904	0.0249032846153157	0.0225138982726739	0.0181697023167219	0.0136609854949669	0.0121135163669742	0.014408338225588
New_Guinea SD	0.0890975299550152	0.0842301251797607	0.0792178489696074	0.0736641677047341	0.0671594485979282	0.0596008898697251	0.0511674521224338	0.0422240170397645	0.0331812675142621	0.0254299321136062	0.0225457007202054	0.0268237812941903	0.034037094432116	0.0390786810653157	0.0406181672896982	0.0401664674846312	0.0399927675349423	0.0409680658412061
Mesoamerica SD	0.0257631275045289	0.0187178864900728	0.0164812107253495	0.0223842305919749	0.0343001607587918	0.0492203013799705	0.0648668376274174	0.0787357939871162	0.0882653199553832	0.0913017452354271	0.0865857410678889	0.0743419528904328	0.0572427399139474	0.0406792417457969	0.030609493113632	0.0315030760258884	0.0386540722319049	0.0479903934364773
N_Lowland_SA SD	0.111921322285171	0.107282514497114	0.102185138216149	0.0958207981245285	0.0873593968043197	0.0763936956444485	0.0630414035072374	0.0480814037523411	0.0329504357663719	0.0219371014627934	0.0252549283205093	0.0397679138235974	0.0556648469342838	0.0682249344844884	0.0766047903929501	0.0827208213023969	0.0895838817981019	0.0974704711325054
NW_Lowland_SA SD	0.150541353866891	0.139669442413694	0.127930041006379	0.114205838438166	0.097868424572761	0.0797876044778232	0.0631128649695468	0.0525350681168884	0.0499851737421644	0.0520614502101073	0.0541664339731254	0.054016322802996	0.0513832224992821	0.0467970105936129	0.0411427020149635	0.0364198792817815	0.035220746530383	0.0381254553823437
Sava_W_India SD	0.00620248739705606	0.00409602037065882	0.00439841375218349	0.00750428764031471	0.012286775530684	0.0183044382262446	0.0252698645356455	0.0327879406646454	0.0403222428411792	0.0470625286607487	0.0518738127485427	0.0537645121764797	0.052235535975566	0.0476849083553341	0.0415125953320544	0.0349790611537332	0.029039548251698	0.0251114044895353
S_India SD	0.0085115048531172	0.00892906518867892	0.00967136329810452	0.0109786376784981	0.0132686066270033	0.0165765384328679	0.0204467930723038	0.0238711849145664	0.0256672972233344	0.0255092175722659	0.0237276529091089	0.0208724202116767	0.0181457061716282	0.0170301004732513	0.0168873351736634	0.0164690070990293	0.0154340912923534	0.0142507144708916
Ganges_E_Indi SD	0.00996215193394363	0.00691767159299303	0.00791433587724963	0.0127874855691734	0.0199162980428461	0.0287137566212716	0.0388527650554929	0.0499222314748722	0.0614506116687209	0.0728233254374634	0.0831495621605526	0.0908947825786612	0.0944845150273855	0.0928545508935936	0.0864932764809146	0.0768625361267711	0.065847656098237	0.0558662960341165
Chinese_loess SD	0.00723695180791585	0.00653465042035558	0.0067229060136461	0.00783510017868662	0.00968749896229179	0.0121031590893252	0.015015434915449	0.0186404465048965	0.0229918308268338	0.0276609710250012	0.0323906245272664	0.0365700364119769	0.0388300984431191	0.037141136645161	0.0311556397194864	0.0229315384793928	0.0156852378394716	0.0147092261255239
Japanese SD	0.184696194095272	0.167622323188328	0.150274927059077	0.131797262192788	0.11137082376895	0.0892794740360921	0.0676920587877228	0.050106916646986	0.0386317289017512	0.0331345985067326	0.0322763726395578	0.0340554941652883	0.0358963772505133	0.0361379405079019	0.035343175991101	0.0349568839559365	0.036562520951131	0.0402092257222996
Lower-MiddleY SD	0.0109221052239635	0.00879404210265073	0.0082644102269953	0.00985094611053262	0.0133492834979398	0.0182250521629163	0.0238393891749069	0.0293168212659387	0.0336239347179238	0.0358993675826794	0.0357122875762433	0.033336309677034	0.0300883638849522	0.0274478338873639	0.0261092333651014	0.0261920192405151	0.0270668753710159	0.0286050078015116
South trop ch SD	0.0164074485323459	0.013017339693911	0.0101269327626496	0.00883966436288355	0.0113546594405518	0.0176935756740366	0.0263781784905157	0.0362001386983728	0.0458282304601001	0.0540734242720755	0.0596617127431108	0.0616233193553812	0.0592826758862225	0.0525468261222113	0.0423612740395727	0.0308335216922742	0.0220855802940787	0.0229498120978104
NA SD	0.0809467917305336	0.0685832123328493	0.0557119272621779	0.0416794929850676	0.0273556018582627	0.0213125874557977	0.0336477364276112	0.0522922853903463	0.0688255532152485	0.0796532138108896	0.0829295203707789	0.0782121775868238	0.0673634830545416	0.0546060927712949	0.0437646882775398	0.0377464374007277	0.0370323136420871	0.0390839743471241
Southwes amaz SD	0.125774207173705	0.115361452592516	0.1042880086271	0.0916253149336918	0.0769919993552501	0.0614840306437265	0.0481765788042826	0.0417025596083742	0.0432490383179266	0.0477994697641187	0.0500604213341766	0.0477221942386258	0.0415134022061204	0.0341998323542236	0.0287859842127018	0.0282643985012614	0.0333826815416954	0.0417094659365358
C/S_Andes SD	0.0727773302277296	0.06676048268613	0.0609678447564098	0.055024169996668	0.0491091885343267	0.0446462966444592	0.0438008333720804	0.0470055485124615	0.0524694588100196	0.0582765439649155	0.0633094382982524	0.0669266582058461	0.068773937206034	0.0681584149345244	0.0650020480377401	0.0612611199316851	0.059915336729057	0.0611532400816508

par(mfrow=c(4,5), mar=c(0,0,0,0))
for(j in 1:20){
plot(means_matrix[1,2:19], col=adjustcolor("cornflowerblue", alpha=0.8), pch=names(means_matrix[1,2:19]), type="n", xlab="year", ylab="Density", xaxt="n", ylim=c(0,1))
originI <- Pop[per.origin[[j]][, 1], ]
for(h in 1:18){
points(23- c(jitter(rep(4, length(originI[,h])), 4.5)  + h), originI[,h], pch=19, cex=.1)
}
  
lines(means_matrix[j,2:19], col=adjustcolor("cornflowerblue", alpha=0.8), pch=names(means_matrix[j,2:19]), type="b", xlab="year", ylab="Density", xaxt="n")
}
axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))

par(mfrow=c(4,5), mar=c(0,0,0,0))
for_3 <- diff_df(1)
global.means <- for_3[[1]]
global.SD <- for_3[[2]]
for(i in 1:20){
  
plot(4:21, global.means[[i]], col=adjustcolor("cornflowerblue", alpha=0.8),  xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(0,1), xlim=c(0,22), xaxt="n", yaxt="n")
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x= 22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))  
  
polygon(c(4:21,21:4) -3 , c(global.means[[i]] + abs(global.SD[[i]]), rev(global.means[[i]] - abs(global.SD[[i]]))), col="cornflowerblue")
lines(global.means[[i]])
#axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
mtext(names(global.means)[i], 3, line=-1, cex=.5, adj=.3)
}

GAM model using one degree of freedom

par(mfrow=c(4,5), mar=c(0,0,0,0))
for_3 <- diff_df(2)
global.means <- for_3[[1]]
global.SD <- for_3[[2]]
for(i in 1:20){
plot(4:21, global.means[[i]], col=adjustcolor("cornflowerblue", alpha=0.8),  xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(0,1), xlim=c(0,22), xaxt="n", yaxt="n")
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x= 22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
polygon(c(4:21,21:4) -3 , c(global.means[[i]] + abs(global.SD[[i]]), rev(global.means[[i]] - abs(global.SD[[i]]))), col="cornflowerblue")
lines(global.means[[i]])
#axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
mtext(names(global.means)[i], 3, line=-1, cex=.5, adj=.3)
}

GAM model using two degree of freedom

par(mfrow=c(4,5), mar=c(0,0,0,0))
for_3 <- diff_df(3)
global.means <- for_3[[1]]
global.SD <- for_3[[2]]
for(i in 1:20){
plot(4:21, global.means[[i]], col=adjustcolor("cornflowerblue", alpha=0.8),  xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(0,1), xlim=c(0,22), xaxt="n", yaxt="n")
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x= 22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
polygon(c(4:21,21:4) -3 , c(global.means[[i]] + abs(global.SD[[i]]), rev(global.means[[i]] - abs(global.SD[[i]]))), col="cornflowerblue")
lines(global.means[[i]])
#axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
mtext(names(global.means)[i], 3, line=-1, cex=.5, adj=.3)
}

GAM model using three degree of freedom

par(mfrow=c(4,5), mar=c(0,0,0,0))
for_3 <- diff_df(4)
global.means <- for_3[[1]]
global.SD <- for_3[[2]]
for(i in 1:20){
plot(4:21, global.means[[i]], col=adjustcolor("cornflowerblue", alpha=0.8),  xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(0,1), xlim=c(0,22), xaxt="n", yaxt="n")
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x= 22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
polygon(c(4:21,21:4) -3 , c(global.means[[i]] + abs(global.SD[[i]]), rev(global.means[[i]] - abs(global.SD[[i]]))), col="cornflowerblue")
lines(global.means[[i]])
#axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
mtext(names(global.means)[i], 3, line=-1, cex=.5, adj=.3)
}

GAM model using four degree of freedom

par(mfrow=c(4,5), mar=c(0,0,0,0))
for_3 <- diff_df(5)
global.means <- for_3[[1]]
global.SD <- for_3[[2]]
for(i in 1:20){
plot(4:21, global.means[[i]], col=adjustcolor("cornflowerblue", alpha=0.8),  xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(0,1), xlim=c(0,22), xaxt="n", yaxt="n")
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x= 22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
polygon(c(4:21,21:4) -3 , c(global.means[[i]] + abs(global.SD[[i]]), rev(global.means[[i]] - abs(global.SD[[i]]))), col="cornflowerblue")
lines(global.means[[i]])
#axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
mtext(names(global.means)[i], 3, line=-1, cex=.5, adj=.3)
}

GAM model using 5 degree of freedom

par(mfrow=c(4,5), mar=c(0,0,0,0))
for_3 <- diff_df(10)
global.means <- for_3[[1]]
global.SD <- for_3[[2]]
for(i in 1:20){
plot(4:21, global.means[[i]], col=adjustcolor("cornflowerblue", alpha=0.8),  xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(0,1), xlim=c(0,22), xaxt="n", yaxt="n")
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x= 22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))  
  
polygon(c(4:21,21:4) -3 , c(global.means[[i]] + abs(global.SD[[i]]), rev(global.means[[i]] - abs(global.SD[[i]]))), col="cornflowerblue")
lines(global.means[[i]])
#axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
mtext(names(global.means)[i], 3, line=-1, cex=.5, adj=.3)
}

GAM model using 10 degree of freedom

par(mfrow=c(4,5), mar=c(0,0,0,0))
for_3 <- diff_df(12)
global.means <- for_3[[1]]
global.SD <- for_3[[2]]
for(i in 1:20){
plot(4:21, global.means[[i]], col=adjustcolor("cornflowerblue", alpha=0.8),  xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(0,1), xlim=c(0,22), xaxt="n", yaxt="n")
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x= 22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))  
  
polygon(c(4:21,21:4) -3 , c(global.means[[i]] + abs(global.SD[[i]]), rev(global.means[[i]] - abs(global.SD[[i]]))), col="cornflowerblue")
lines(global.means[[i]])
#axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
mtext(names(global.means)[i], 3, line=-1, cex=.5, adj=.3)
}

GAM model using 12 degree of freedom

par(mfrow=c(4,5), mar=c(0,0,0,0))
for_3 <- diff_df(17)
global.means <- for_3[[1]]
global.SD <- for_3[[2]]
for(i in 1:20){
plot(4:21, global.means[[i]], col=adjustcolor("cornflowerblue", alpha=0.8),  xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(0,1), xlim=c(0,22), xaxt="n", yaxt="n")
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x= 22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))   
  
polygon(c(4:21,21:4) -3 , c(global.means[[i]] + abs(global.SD[[i]]), rev(global.means[[i]] - abs(global.SD[[i]]))), col="cornflowerblue")
lines(global.means[[i]])
#axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
mtext(names(global.means)[i], 3, line=-1, cex=.5, adj=.3)
}

GAM model using 17 degree of freedom

Productivity

# Load patricks productivity PCA data
load('Productivity_ALL.RDATA')
# Load origin shapefiles
origins <- readShapePoly('Origins_updated.shp')
origin.time.region <- c(2, 2, 1, 1, 1, 2, 2, 1, 2, 2, 
                        2, 2, 1, 2, 2, 2, 2, 2, 2, 2) # 1 = early; 2 = middle
# Extract the data
prod.origin <- extract(Productivity.ALL, origins)

NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files

# Mean and SD per region
means <- lapply(prod.origin, colMeans, na.rm = TRUE)
sds <- lapply(prod.origin, sd, na.rm = TRUE)
names(means) <- origins@data$CONTINENT
ymax <- max(unlist(means))
ymin <- min(unlist(means))
time <- 4:21
# Plot
#pdf("productivity.pdf", 20, 30) 
par(mfrow = c(5, 4), mar = c(2, 2, 2, 0))
for (i in 1:length(means)) {
  plot(y = means[[i]], x = time, xlim = c(21, 4), ylim = c(ymin, ymax),
       main = names(means)[i], cex.main = 1, cex.lab = 1, cex.axis = 1,
       ylab = "Productivity (PCA axis)", xlab = "Thousand of years ago (k)",
       pch = 20, lwd = 1, type = "l", 
       col = c("purple", "green")[origin.time.region[i]])
  up <- sds[[i]] + means[[i]]
  down <-  means[[i]] - sds[[i]]
  lines(up ~ time, lty = 2)
  lines(down ~ time, lty = 2)
  
}

#dev.off()

Compare rates between origins and not-origins

Compare rates between different time periods

Setup final figure

Frame in the layout

a <- layout(matrix(c(
    1, 1, 1, 1, 1, 1, 1, 1,
    3,  6, 7, 8, 9, 10, 11, 4, 
    3,  5, 5, 5, 5, 5, 5,   4, 
    3,  12, 13, 14, 15, 16, 17, 4,
    2, 2, 2, 2, 2, 2, 2, 2
    ), 5, 8, byrow=TRUE), width=c(1, 1, 1, 1, 1, 1, 1, 1), height=c(0.5, 1, 1.5, 1, 0.5))
layout.show(a)

Make blank template plots

frameplot <- function(){
    plot(21:0,rep(0, 22), xlim=c(17,4), ylim=c(0, 2.25), type="n", xaxt="n", yaxt="n", xlab="", ylab="")
}
frameplot_bottom <- function(){
    plot(21:0,rep(0, 22), xlim=c(17,4), ylim=c(-0.25, 2), type="n", xaxt="n", yaxt="n", xlab="", ylab="")
}

frameplot()

frameplot_bottom()

map for final figure

Make the map for the center panel (#5 on layout panel)

d <- readPNG("earth.png")

png(file=paste("40962.png",sep=""),width=2000,height=1000, bg="transparent")
par(mar=c(0,0,0,0))
plot(seq(-180, 180, length.out = 19), seq(-90, 90, length.out = 19), type="n",xlim=c(-180, 180),ylim=c(-90, 90), xaxt="n")
rasterImage(d, -180, -90, 180, 90, interpolate=TRUE, col=d)
polygon(x=c(-180,-180, 180,180), y=c(-90, 90, 90, -90), col=adjustcolor("white", alpha=0.1))
#rasterImage(d, -13.5, -13.5, 375, 375, interpolate=TRUE, col=d)
plot(origins_subset, add=TRUE, col=adjustcolor("white", alpha=.8), xaxt="n", border="white", lwd=4) #still need to reproject!!!
dev.off()

null device 
          1 

Trend through time panel

Setup the plot template for small panel plots (#6-17 on layout panel)

###################
type_number <- 8
complex_figure <- function(type_number, i, means, sds){
                        
if(i < 6)   polygon(x=c(12,12,8.2,8.2), y=c(-1,3,3,-1), col=adjustcolor("cornflowerblue", alpha=0.4), border=NA)                    
if(i > 5)   polygon(x=c(8.2,8.2,4.2,4.2), y=c(-1,3,3,-1), col=adjustcolor("limegreen", alpha=0.4), border=NA)
                                    
    match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 9]
    maxer <- max(match, na.rm=TRUE)
    j <- ecdf(maxer-match)
    print(j)
    
x_seq <- rev(c(0,seq(0, maxer, length.out=100)))
y_seq <- -c(0, j(seq(0, maxer, length.out=100)))
#lines(x_seq, y_seq, type="l", ylim=c(-1,1))
#polygon(c(0, x_seq), c(0, y_seq), border="black", col=adjustcolor("cornflowerblue", alpha=0.5))
#abline(v= maxer - quantile(j)[2])
    
    break_one_1 <- maxer
            break_two_1 <- maxer - quantile(j)[2]
                
#   polygon(x=c(break_two_1, break_two_1, break_one_1, break_one_1), y=c(0, 1, 1, 0), col=adjustcolor("cornflowerblue", alpha=0.5), border=NA)
            
    match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 10]
    maxer <- max(match, na.rm=TRUE)
    j <- ecdf(maxer-match)
    print(j)
    
x_seq <- rev(c(0,seq(0, maxer, length.out=100)))
y_seq <- 2+c(0, j(seq(0, maxer, length.out=100)))
#lines(x_seq, y_seq)
#polygon(c(0, x_seq), c(2, y_seq), border="black", col=adjustcolor("limegreen", alpha=0.5))
    
    break_one_2 <- maxer
            break_two_2 <- maxer - quantile(j)[2]
                
#   polygon(x=c(break_two_2, break_two_2, break_one_2, break_one_2), y=c(1, 2, 2, 1), col=adjustcolor("limegreen", alpha=0.5), border=NA)
            
    
        #abline(v=11)
    type <- 1
        
        if(type == 1){
    x <- c(means[[i]] , means[[i]]  + abs(sds[[i]]), means[[i]]  - abs(sds[[i]]))
    scaled <- scale(x , center=FALSE)
    meanss <- scaled[1:18]
    sdss_plus <- scaled[19:36]
    sdss_minus <- scaled[37:54]
    #abline(v=10, col="red")
    length(scaled)
    #lines(4:21, means[[i]] + sds[[i]])
    #polygon(x=c(4:21, 21:4), y=c(sdss_plus, rev(sdss_minus)), col=adjustcolor("firebrick", alpha=1), border="white")
    polygon(x=c(21:4,4:21), y=c(sdss_plus, rev(sdss_minus)), col=adjustcolor("firebrick", alpha=1), border="white") 
    }
    
    if(type == 2){
    x <- c(means[[i]] , means[[i]]  + abs(sds[[i]]), means[[i]]  - abs(sds[[i]]))
    scaled <- x + 1 #scale(x , center=FALSE)
    meanss <- scaled[1:18]
    sdss_plus <- scaled[19:36]
    sdss_minus <- scaled[37:54]
    #abline(v=10, col="red")
    length(scaled)
    #lines(4:21, means[[i]] + sds[[i]])
    polygon(x=c(4:21, 21:4), y=c(sdss_plus, rev(sdss_minus)), col=adjustcolor("firebrick", alpha=1), border="white")    
    }
if(type == 3){
    x <- c(means[[i]] , means[[i]]  + abs(sds[[i]]), means[[i]]  - abs(sds[[i]]))
    scaled <- x #scale(x , center=FALSE)
    meanss <- scaled[1:18]
    sdss_plus <- scaled[19:36]
    sdss_minus <- scaled[37:54]
    #abline(v=10, col="red")
    length(scaled)
    #lines(4:21, means[[i]] + sds[[i]])
    polygon(x=c(4:21, 21:4), y=c(sdss_plus, rev(sdss_minus)), col=adjustcolor("firebrick", alpha=1), border="white")    
    }
    
    
means_long_y <- c(1,1,1,1,1, meanss)
means_long_x <- c(0:4, 4:21)
 
            break_one <- break_one_2
            break_two <- break_two_2
        #       polygon(x=c(break_one, break_one, 22, 22), y=c(1, 2, 2, 1), col=adjustcolor("white", alpha=0.8), border=NA)
        #       polygon(x=c(break_two, break_two, break_one, break_one), y=c(1, 2, 2, 1), col=adjustcolor("white", alpha=0), border=NA)
            #   polygon(x=c(-1,-1, break_two , break_two), y=c(1.9, 3.1, 3.1, 1.9), col=adjustcolor("white", alpha=0.8), border=NA) 
                #abline(v= break_one, col="white")
                #abline(v= break_two, col="white")
                
                break_one <- break_one_1
            break_two <- break_two_1
        #       polygon(x=c(break_one, break_one, 22, 22), y=c(0, 1, 1, 0), col=adjustcolor("white", alpha=0.8), border=NA)
        #       polygon(x=c(break_two, break_two, break_one, break_one), y=c(0, 1, 1, 0), col=adjustcolor("white", alpha=0), border=NA)
            #   polygon(x=c(-1,-1, break_two , break_two), y=c(-1.1, .1, .1, -1.1), col=adjustcolor("white", alpha=0.8), border=NA) 
                #abline(v= break_one, col="white")
                #abline(v= break_two, col="white")
                
#lines(x=c(break_one_2, break_one_2), y=c(1,3), col="white")
#lines(x=c(break_one_1, break_one_1), y=c(1,-1), col="white")
#lines(x=c(break_two_2, break_two_2), y=c(1,3), col="white")
#lines(x=c(break_two_1, break_two_1), y=c(1,-1), col="white") 
#lines(4:21, meanss)
    lines(21:4, meanss)
    
}

frameplot()
complex_figure(7, 1, global.means, global.SD) 

Empirical CDF 
Call: ecdf(maxer - match)
 x[1:7] =      0,    0.5,      1,  ...,      5,    5.5
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:4] =      0,    3.5,      4,    4.5

axis(1)
axis(2)

Assemble the figure

Assemble the figure

quartz(width=8, height=8)
layout(matrix(c(
    1, 1, 1, 1, 1, 1, 1, 1,
    3,  6, 7, 8, 9, 10, 11, 4, 
    3,  5, 5, 5, 5, 5, 5,   4, 
    3,  12, 13, 14, 15, 16, 17, 4,
    2, 2, 2, 2, 2, 2, 2, 2
    ), 5, 8, byrow=TRUE), width=c(1, 1, 1, 1, 1, 1, 1, 1), height=c(0.5, 1, 1.5, 1, 0.5))
par(mar=c(0,0,0,0))
# 1-4 label margins
blankplot <- function(){
    
    plot(0,0, xlim=c(4,21), ylim=c(1, 1.25), bty="n", type="n", xaxt="n", yaxt="n", xlab="", ylab="")
}
blankplot()
blankplot()

blankplot()
blankplot()

origins <- readShapePoly('Origins_updated.shp')
origin.time.region <- c(2, 2, 1, 1, 1, 2, 2, 1, 2, 2, 
                        2, 2, 1, 2, 2, 2, 2, 2, 2, 2) # 1 = early; 2 = middle
as.character(origins$CONTINENT)

 [1] "W_African_Sav" "Sudanic_Savan" "West Africa T" "Ethipian plat" NA             
 [6] "E_North_Ameri" "New_Guinea"    "Mesoamerica"   "N_Lowland_SA"  "NW_Lowland_SA"
[11] "Sava_W_India"  "S_India"       "Ganges_E_Indi" "Chinese_loess" "Japanese"     
[16] "Lower-MiddleY" "South trop ch" NA              "Southwes amaz" "C/S_Andes"    

#subset_order <- c(1, 2, 3, 5, 6, 8, 9, 10, 11, 12, 17, 18)
subset_order <- c(8, 10, 9, 5, 18, 7, 6, 20, 1, 2, 13, 16)
origins_subset <- origins[subset_order,]
origins_subset$CONTINENT

 [1] Mesoamerica   NW_Lowland_SA N_Lowland_SA  <NA>          <NA>          New_Guinea   
 [7] E_North_Ameri C/S_Andes     W_African_Sav Sudanic_Savan Ganges_E_Indi Lower-MiddleY
18 Levels: C/S_Andes Chinese_loess E_North_Ameri Ethipian plat Ganges_E_Indi ... West Africa T

d <- readPNG("earth.png")
png(file=paste("40962.png",sep=""),width=2000,height=1000, bg="transparent")
par(mar=c(0,0,0,0))
plot(seq(-180, 180, length.out = 19), seq(-90, 90, length.out = 19), type="n",xlim=c(-180, 180),ylim=c(-90, 90), xaxt="n")

rasterImage(d, -180, -90, 180, 90, interpolate=TRUE, col=d)
polygon(x=c(-180,-180, 180,180), y=c(-90, 90, 90, -90), col=adjustcolor("white", alpha=0.1))
#rasterImage(d, -13.5, -13.5, 375, 375, interpolate=TRUE, col=d)
plot(origins_subset, add=TRUE, col=adjustcolor("white", alpha=.8), xaxt="n", border="white", lwd=4) #still need to reproject!!!
dev.off()

quartz 
     2 

d <- readPNG("40962.png")
dim(d)

[1] 1000 2000    4

par(mar=c(0,0,0,0))

plot(0:360,0:360,type="n",xlim=c(20,360),ylim=c(65,295), yaxt="n", xaxt="n")
rasterImage(d, -28.5, -13.5, 388, 375, interpolate=TRUE, col=d)

axis(2, label=seq(-90, 90, length.out = 19), at=seq(1, 360, length.out = 19), las=1)
mtext("latitude", 2, line=4, at=180)

abline(h=seq(1, 360, length.out = 19), col=adjustcolor("grey10", alpha= 0.4), lwd=1)
abline(h=180, col=adjustcolor("white", alpha= .5), lwd=1)

load('PopD_all_December.rdata')
# Extract the data
prod.origin <- extract(PopD.ALL, origins_subset)

NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files

library(matrixStats)

matrixStats v0.51.0 (2016-10-08) successfully loaded. See ?matrixStats for help.

# Mean and SD per region
means <- lapply(prod.origin, colMeans, na.rm = TRUE)
sds <- lapply(prod.origin, colSds, na.rm = TRUE)
## new values from Bruno's GAM model (produced in script called Population_Trend_per_y.R)
means <- global.means
sds <- global.gams
names(means) <- origins_subset@data$CONTINENT
ymax <- max(unlist(means))
ymin <- min(unlist(means))
time <- 4:21
#plot(origins)
#means[[1]] +
#sds[[1]]
#scale(as.numeric(means[[1]]), center=FALSE)
name_vector <- as.character(origins_subset@data$CONTINENT)
for(i in 1:12){
    
    if(i > 6){frameplot()}else{frameplot_bottom()}
    
        ## customize polygons for each graph
    if(i == 1){ #mesoamerica  #values from Larson
        
            complex_figure(3, i, means, sds)
                
    
        }
    
    
    #########
    if(i == 2 ){ #NW lowlands SA  #values from Larson
        
        complex_figure(6, i, means, sds)
    
        }
        
        #########
    if( i == 3){ #NW lowlands SA  #values from Larson
        
        complex_figure(6, i, means, sds)
        
        }
        #########
    if(i == 4){ #Fertile crescent aka Southwest asia  #values from Larson
        
        
    complex_figure(8, i, means, sds)
                
        }
        
        #########
    if(i == 5){ #loess plateau  #values from Larson
        
        complex_figure(2, i, means, sds)
            
        }
        
        
        #########
    if(i == 6){ #new guinea  #values from Larson
        
        complex_figure(4, i, means, sds)
        
        }
#########
    if(i == 7){ #Eastern N.A.  #values from Larson
        
        complex_figure(5, i, means, sds)
        
            }
        #########
    if(i == 8){ #Andes  #values from Larson
        
        complex_figure(6, i, means, sds)
        
                }
#########
    if(i == 9){ #W. African Sav  #values from Larson
        
        complex_figure(1, i, means, sds)
        
            }
#########
    if(i == 10){ #Sudanic sav  #values from Larson
        
        complex_figure(1, i, means, sds)
        
                }
#########
    if(i == 11){ #Ganges  #values from Larson
        
        
        complex_figure(7, i, means, sds) 
        
        }
#########
    if(i == 12){ #loess  #values from Larson
        
        complex_figure(2, i, means, sds)
         
                }
        
        
        #lines(4:21, means[[i]])
        
        abline(h = 1, col=adjustcolor("forestgreen", alpha=.5), lty=2)
        
    # add axes to some locations
    if(i == 1 | i == 7){axis(2, at=seq(0,2, by=0.25), label=seq(0,2, by=0.25), las=1)}
    if(i == 6 | i == 12){axis(4, at=seq(0,2, by=0.25), label=seq(0,2, by=0.25), las=1)}
    #if(i == 6 | i == 12){axis(4, at=seq(2,3, by=0.25), label=seq(0,1, by=0.25), las=1)
    #   axis(4, at=seq(-1,0, by=0.25), label=rev(seq(0,1, by=0.25)), las=1)
    #   }
    if(i > 6){axis(1)} else{axis(3)}
    
    # add text 
    if(i < 7){polygon(x=c(-30, -30, 30, 30), y=c(-0.1, -0.5, -0.5, -0.1), col="black")
    mtext(name_vector[i], 1, line=-1.2, col="white", cex=0.5)}
    
    if(i > 6){polygon(x=c(-30, -30, 30, 30), y=c(2.1, 2.5, 2.5, 2.1), col="black")
    mtext(name_vector[i], 3, line=-1.2, col="white", cex=0.5)}
    
    # add axis labels
    if(i == 1 | i ==  7){mtext("scaled density potential", 2, line=4, at=1)}
    if(i ==  3){mtext("Thousand years before present", 3, line=3.5, at =5)}
    if(i ==  9){mtext("Thousand years before present", 1, line=3.5, at =5)
        
        }
    
}

Empirical CDF 
Call: ecdf(maxer - match)
 x[1:4] =      0,      4,      7,      8
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:5] =      0,      1,      4,      7,      8
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:8] =      0,      1,      2,  ...,      5,      6
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:9] =      0,      1,      2,  ...,      6,    6.5
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:8] =      0,      1,      2,  ...,      5,      6
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:9] =      0,      1,      2,  ...,      6,    6.5
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:6] =      0,    0.5,      1,  ...,      2,      8
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:8] =      0,    0.5,   0.75,  ...,      5,      7
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:11] =      0,    0.5,      2,  ...,   7.75,    9.5
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:9] =      0,    0.5,      1,  ...,    4.5,    7.5
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:1] =      0
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:2] =      0,      4
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:3] =      0,    0.5,   1.75
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:3] =      0,      1,   1.25
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:8] =      0,      1,      2,  ...,      5,      6
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:9] =      0,      1,      2,  ...,      6,    6.5
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:11] =      0,   0.25,   1.25,  ...,   6.75,   7.75
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:11] =      0,   2.25,   2.75,  ...,   5.75,   6.25
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:11] =      0,   0.25,   1.25,  ...,   6.75,   7.75
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:11] =      0,   2.25,   2.75,  ...,   5.75,   6.25
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:7] =      0,    0.5,      1,  ...,      5,    5.5
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:4] =      0,    3.5,      4,    4.5
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:11] =      0,    0.5,      2,  ...,   7.75,    9.5
Empirical CDF 
Call: ecdf(maxer - match)
 x[1:9] =      0,    0.5,      1,  ...,    4.5,    7.5

saveToPDF <- function(...) {
    d = dev.copy(pdf,...)
    dev.off(d)
}
saveToPNG <- function(...) {
    d = dev.copy(png,...)
    dev.off(d)
}
## Try them out
saveToPDF("my.pdf", height=8,width=8)

quartz 
     2 

saveToPNG("my.png", height=8, width=8, units="in", res=300)

quartz 
     2 

dev.off()

null device 
          1